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Abstract 

The BGK model kinetic equation is applied to spatially inhomogeneous states 
near steady uniform shear flow. The shear rate of the reference steady state 
can be large so the states considered include those very far from equilibrium. 
The single particle distribution function is calculated exactly to first order in 
the deviations of the hydrodynamic field gradients from their values in the 
reference state. The corresponding non-linear hydrodynamic equations are 
obtained and the set of transport coefficients are identified as explicit functions 
of the shear rate. The spectrum of the linear hydrodynamic equations is 
studied in detail and qualitative differences from the spectrum for equilibrium 
fluctuations are discussed. Conditions for instabilities at long wavelengths are 
identified and discussed. 

PACS numbers: 05.20Dd, 51.10+y, 83.50Ax, 83.20-d 
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I. INTRODUCTION 



Nonequilibrium phenomena are well understood for states near equilibrium. However, 
even the qualitative features of transport and fluctuations far from equilibrium are poorly 
understood, due to the complexity of the physical states and the lack of an adequate con- 
trolled theoretical description in general. It is useful in this context to consider more re- 
strictive conditions and specialized states for which greater progress in this difficult area 
can be made. Here we consider a low density simple atomic gas for which transport prop- 
erties are described by the non-linear Boltzmann equation. For states near equilibrium, the 
Chapman-Enskog expansion of the distribution function about the local equilibrium distri- 
bution in terms of gradients of the hydrodynamic fields provides approximations to a normal 
solution (one for which all space and time dependence occurs through the hydrodynamic 
variables) From this solution the corresponding hydrodynamic equations are obtained, 
including explicit expressions for the associated transport coefficients. In principle, this 
method applies to states far from equilibrium as well, although calculation of the Chapman- 
Enskog expansion to higher orders in the gradients is prohibitively difficult and questions 
of convergence remain unresolved. An alternative approach is to expand in small gradients 
about a more relevant reference state than local equilibrium. For example, consider states 
near a homogeneous reference steady state. Deviations of the hydrodynamic variables from 
their values in this state are characterized by small relative spatial gradients. A modified 
Chapman-Enskog expansion can be implemented to obtain the distribution function and 
hydrodynamic equations to leading order in these gradients. Since the reference state can 
be far from equilibrium the form of the hydrodynamic equations and the dependence of the 
transport coefficients on parameters of the steady state will be quite different from those for 
states near equilibrium. 

In practice this program has not been carried out since determination of a non-trivial 
reference steady state far from the equilibrium from the Boltzmann equation is exceptionally 
rare @. In contrast, exact results have been obtained for a number of physically interesting 
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special steady states using a model kinetic equation designed to preserve the essential fea- 
tures of the Boltzmann equation while admitting more practical analysis 0]. In some cases 
these results have been compared with those from Monte Carlo simulation far from equilib- 
rium, indicating that the kinetic model provides both a qualitative and semi- quantitative 
representation of the underlying Boltzmann equation |5||. The objective here is to use this 
kinetic model approach for the special class of states at or near uniform shear flow. Uniform 
shear flow is a particularly well-studied nonequilibrium state using numerical methods (both 
Monte Carlo at low density || and molecular dynamics simulation at high density |||7|]), il- 
lustrating rheological properties normally associated with complex molecular fluids [§ . More 
importantly, it is one of the special cases for which an exact solution to the model kinetic 
equation has been obtained [p|,|T0|], providing the necessary reference state for the modified 
Chapman-Enskog expansion developed here. We look for solutions to the kinetic equation as 
an expansion in small spatial gradients relative to an exact solution for local uniform shear 
rather than around local equilibrium. The heat and momentum fluxes are calculated from 
this solution to linear order in these gradients, and a closed set of generalized hydrodynamic 
equations is obtained. This is analogous to the non-linear Navier-Stokes equations near 
equilibrium, except that here the reference state is non-Maxwellian and a complex function 
of the shear rate. The associated nonlinear transport coefficients are identified and examples 
calculated for arbitrary values of the shear rate. 

To expose the physical content of these hydrodynamic equations and their differences 
from those for states near equilibrium, the hydrodynamic modes are calculated from the 
associated linearized equations. In addition to their dependence on the shear rate, these 
modes have a more complex wavevector dependence than those from the Navier-Stokes 
equations due to the broken symmetry of the reference state. More surprisingly, a new long 
wavelength instability is found such that the hydrodynamic modes are growing initially for 
any finite shear rate and sufficiently long wavelength [[LI]]. A more detailed study of this 
instability and its verification via computer simulation is described elsewhere JT2] . 
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II. KINETIC THEORY AND UNIFORM SHEAR FLOW 



In this section the kinetic theory is defined, the associated macroscopic conservation laws 
obtained, and the special solution for steady uniform shear flow is described. In the next 
section, this solution is generalized to a local reference state for an expansion to describe a 
class of states near uniform shear flow and to obtain the associated hydrodynamic equations. 

Exact or even approximate solutions to the Boltzmann equation far from equilibrium 
are exceptionally rare, due to the complexity of the nonlinear collision operator. Therefore 
kinetic models have been introduced to replace the Boltzmann collision operator with a 
simpler, more tractable operator. The best studied of these is a single relaxation time model 
due to Bhatnager, Gross and Krook (BGK model) Q. The essential qualities of this model 
are its preservation of the exact equilibrium solution and all five conservation laws. The 
BGK-Boltzmann kinetic equation is given by 

{e + V ' Vr ) f + m ~ 1Vv ' (Fea;</) = ~" (/ " (2,1) 

where F ext is an external force. The parameter v in the fl2.1|) is a collision frequency which 
depends on an interaction law. This frequency is a function of the density and temperature. 
At low density, it can be written as 

z/~n(r,t)T 6 (r,t), (2.2) 

when the potential has a form V(r) ~ r~', with 6 = 1/2 — 2/1. In the case of Maxwell 
molecules (I = 4), b is zero so v becomes independent of temperature. For the hard sphere 
case / — > oo resulting in a value for b of 1/2. Finally, fg(r,v,t) is the local equilibrium 
distribution 

£(r,v,f) = n(r,t)(^^) 3 / 2 exp (-i/?(r,t)m(v - u(r,t)) 2 ) , (2.3) 

where n(r, i), T(r,t) = [k B [3(r, t)] _1 , and u(r, t) are the density, temperature, and flow 
velocity of the nonequilibrium state. These hydrodynamic fields are defined such that 
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dv 



v 

,2 



(f(r,v,t)-f e (r,v,t)) = 0, 



(2.4) 



which assures that the BGK equation yields the correct conservation laws and equilibrium 
stationary state in the absence of driving forces. More explicitly ( |2.4j ) gives 



n(r,t) = J dvf(r,v,t), n(r, t)u(r, t) = J dv vf(r, v, t), 

3 /" 1 

-n(r,t)k B T(r,t) = J dv-mc 2 f(r,v,t), 



(2.5) 



where c = v — u. 

The macroscopic conservation laws are obtained by taking moments of Q2.1|) and using 
the definitions (|2.5| ) 



D t n(r,t) + n(r,t)V- u(r,t) = 0, 



(2.6) 



D t T{r,t) + 



3fcsn(r, £) 



[V ■ q(r, t) + %(r, t)^(r, t)] = s(r, t) 



D t Ui(r, t) + [mn(r, t)] 1 = 0, 



(2.7) 



(2.8) 



where _D t = 9t + u ■ V is the material derivative. The heat and momentum fluxes, q(r, t) 
and tij(r, t), are linear functionals of /(r, v, t) given by 



q(r,t) = J dv^mc 2 cf(r,v,t), Uj{v,t) = J dv m c i c i /(r, v, £). 



(2.9) 



The inhomogeneous term on the right side of the temperature equation, s(r, t), is due to the 
external force F exf , introduced to serve as a thermostat. Several thermostats that have been 
used in both theory and computer simulations. We choose here a force that is proportional 
to the relative velocity c = v — u(r, t) 



F«rt(r,c,t) = -m\(n(r,t),T(T,t))c. 



(2.10) 



The resulting source term s(r, t) in the equation for the temperature becomes 

s(r,t) = -2T(r,t)X(n(r,t),T(r,t)). (2.11) 

The proportional "constant" X(n(r, t), T(r, t)) is determined by requiring stationarity of the 
system in the uniform shear flow state (see below), and may depend on the local density, 
temperature, and shear rate. In Appendix ^ a different thermostat is considered for com- 
parison. The primary changes are the degree to which the external force compensates for 
viscous heating away from the state of uniform shear flow. 

The fact that the fields are functionals of /(r, v, t) makes the BGK-Boltzmann equation 
highly non-linear and difficult to solve in general. However, in many cases an implicit solution 
can be given as an explicit function of the velocity and functional of the fields. Then, the 
fields must be determined self-consistently from the above macroscopic conservation laws. 
One of the cases for which an exact solution is known is uniform shear flow [P,p|,|10[. The 
uniform shear state is a planar flow whose x-component of the flow velocity has a gradient 
along the where a is a constant shear rate. In addition, 

the density n s , temperature T s , heat flux, and momentum flux are spatially constant. This 
state is generated by a periodic boundary condition in the local Lagrangian frame. The 
viscous heating induced by these boundary conditions is compensated by the external force. 
It is easily verified that this macroscopic state is an exact stationary solution to the above 
conservation laws ( |2.6|) - (|2.8| ) if A(n(r, t), T(r, t)) is chosen to be 



X s ^X(n s ,T s ) = -p^. (2.12) 
6n s k B ls 

Due to the simplicity of uniform shear flow state at the macroscopic level, it has been studied 



extensively as a prototype of nonequilibrium states far from equilibrium in theory [|10[|14| -|T7 
and in computer simulations []6|,|7|, p~8| ? p~9|] . To obtain the distribution function for this state it 
is useful to express the kinetic equation ( |2.1| ) in terms of the velocity in the local rest frame, 
defined by v[ = Vi — a^rj. In this frame, the flow field vanishes and the macroscopic state 
becomes spatially homogeneous. Consequently, we look for a stationary solution to Eq. Q2.1p 
of the form /(r, v, t) = / s (v') 



d 

L(v , a)f s (v ) = -v s f si >(v), L(v , a) = a^— + A s v -V v > + 3A S - u s . (2.13) 

The subscript s denotes the stationary state value and f S £ is the corresponding local equilib- 
rium distribution function with the hydrodynamic fields for uniform shear flow. The solution 
to Eq. flnSD is 



fs(-v')=v s J o dre TL f e (V) = u s J o dte-^-^Ue^A^-tWj). (2.14) 

The second equality follows from the property for an arbitrary function, X(v) 

e tL X(V) = e( 3A -^X(e A *%(-t)^), K l3 (t) = 6 tJ - a tJ t. (2.15) 

where use has been made of exp(atv x d Vx )X (v x ) = X(e at v x ) and exp(atv y d Vx )X (v x ) = X{v x + 
atVy). To determine X s as an explicit function of a, the component of the momentum flux 
txy,s{o) can be calculated from Eqs. Q2.9p and Q2.14Q to get the self-consistent equation 



3A S (2A S + z/ s ) 2 = u s a 2 . (2.16) 

This has one real solution and two complex conjugate solutions. The physically relevant 
real value is 



X 8 (a) = ^sinhtcosh-^l + g^)] 2 . (2.17) 
3 o ui 



With X s known, the velocity distribution given by (|2.14p is completely determined. 

Any transport property of interest now can be calculated by integration. A detailed dis- 



cussion can be found in [10| and comparison with Monte Carlo simulations of the Boltzmann 



equation for shear flow is given in [HI]. Only the transport properties associated with the 
heat and momentum fluxes are considered further here. These can be calculated directly 
from ( |2.9| ) and ( |2.14 ) with the results 



q(a) = 0, tij,s( a ) = (Ps + ^a 2 ^i(a))^- - r){a){aij + a^) - ^i(a)a ik a jk (2.18) 

Thus, the heat flux vanishes but the momentum flux describes non-trivial rheological effects 
in terms of the hydrostatic pressure, p = uksT, the shear viscosity, 77(a) = — a~ 1 t x .j / (a), and 
the viscometric function, ^i(a) = a~ 2 \t yy (a) — t xx (a)] , where 



V(a) = , 9W ?, „ Ps, *i(a) = ~ 2 , 9 w ( ?L T P- ( 2 - 19 ) 
(2A s (a) + i/ s ) 2 a 2 (2A s (a) + z/ s ) 

In general there is a second independent viscometric function, ^2(0) = a~ 2 [t 2 «(a) — £ y y(a)], 
which vanishes for our kinetic model. The magnitudes of these transport coefficients are 
monotonically decreasing functions of the shear rate, and have been discussed in detail 
elsewhere fTOj . 



III. HYDRODYNAMICS NEAR UNIFORM SHEAR FLOW 

In this section we consider states that deviate from uniform shear flow by small spatial 
gradients. A solution to the BGK-Boltzmann equation (|2.1|) is obtained by a variant of 
the Chapman-Enskog method whereby the distribution function is expanded about a local 
uniform shear flow reference state in terms of the small spatial gradients of the hydrodynamic 
fields relative to those of uniform shear flow. This is analogous to the usual Chapman-Enskog 
expansion about a local equilibrium distribution. The solution obtained in this way can be 
used to calculate the heat and momentum fluxes in terms of the hydrodynamic variables, so 



that Eqs. (|2.6p - (|2.8|) become a closed set of hydrodynamic equations. The analysis here is 
carried out to first order in the gradients. For small shear rate the usual Navier Stokes results 
are recovered, where the heat flux is given by Fourier's law and the momentum flux is given 
by Newton's viscosity law. However, for large shear rates these fluxes and the corresponding 
hydrodynamic equations are more complex. 

To construct the Chapman-Enskog expansion we look for solutions of the form 

f(r,v,t) = f(V,y a (r,t)), (3.1) 

where y a (r, t) are the hydrodynamic fields, and v[ = Vi — a^rj. This representation expresses 
the fact that the space dependence of the reference shear flow is completely absorbed in the 
relative velocity variable, v', and all other space and time dependence occurs entirely through 
a functional dependence on the hydrodynamic variables, y a (r,t). This is an example of a 
"normal" solution, which is expected to result from a wide class of initial conditions at long 
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times and large space scales. It is essential for a hydro dynamic description since the velocity 
average of any quantity becomes a functional of the hydrodynamic fields. Approximate 
solutions to the BGK equation are obtained by expanding ( |3.1| ) in a formal uniformity 
parameter, e, that measures the spatial gradients of the fields y a (r,t) 

/(V, y (r, f)) = / (0) (v\ y a (v, t)) + e/« (V, y a (r, *)) + ■■•■ (3.2) 

This expansion leads to a corresponding expansion for the heat and momentum fluxes when 
substituted into 



ttf=fg ) +etg ) + -" q=q (0) +6q« + ---, (3.3) 

tf) = J dvmCiCjfM, q (r) = J dv^mc 2 cf {r) . (3.4) 

Finally, use of this expansion for the fluxes in the conservation laws ( [2.6|) -( plf ) identifies an 
expansion for the time derivatives of the fields in powers of the uniformity parameter 

a e<°> am 

S = ir + f aT + "" (3 ' 5) 

These results provide the basis for generating the Chapman-Enskog solution to the BGK- 
Boltzmann equation. 

The BGK-Boltzmann equation in terms of the variable v' is obtained directly from (|2.1|) 



(lft + ^ + aijTj ^. ~ L{v '' g) + X6u ' Vv ') /(r ' v '' = z// ^ r ' v '' ^' 



(3.6) 



where 5u = u — u s . Also, the operator L(v',a) is defined by (|2.13| ) except with z/ s ,A s 
replaced by u, A as functions of n(r,t) and T(r,t). The form of A is still to be determined 
at this point. The spatial gradient in the second term on the left side is taken at constant 
v', and consequently this term is of first order in the uniformity parameter. Substituting 
the expansions ( p.2[ ) and ( |3.5| ) into ( |3.6| ) and equating coefficients of each degree in the 



uniformity parameter leads to the equations for and /W 



;)j L(v',a) + \5u-V v ,\fW=vf £ , (3.7) 



(^-L(v',a)+\6u-Vv) / (1) = -( 



To lowest order in the expansion the conservation laws give 



+ W + W)«T / (0) - (3-8) 



-7— = 0, — — h dijOUj = 0, - nk B — — + dijty + 3Xnk B T = 0, (3.9) 

.(0) 



where t\, is defined by (|3.4|). The parameter A of the external force is now chosen to impose 



^ = ie 



X(n{r,t),T{r,t)) = -at^(r,t)/3n{r,t)k B T{r,t). (3.10) 



The solution to ( |3.7| ) is obtained in a way similar to that for ( |2.13| ) with the result 
/"V, faM) = i/(r, i)n(r, i) ^ f * e -*W 

6XP (' 2t B r(r,t) e2A(r '" r|A,j( " T>( "' " fa ' (r ' t))|2 ) ' (3 - U) 
where A^-(i) is defined in ( |2.15| ). The relationship of ( |3.11|) to ( |2.14| ) is analogous to the 
relationship of the local equilibrium distribution to the strict equilibrium distribution, where 
the former is obtained from the latter by replacing the hydrodynamic fields with their actual 
nonequilibrium values. To determine A(r, t) = A(n(r, t),T(r, t)), the momentum flux t^j (r, t) 
is calculated using ( |3.4j ) and ( |3.11|) to obtain 



*g(M) = n(T,t)k B T(r t t) t (3.12) 

(2A(r,t) +v[T,t)) 

where v(r,t) = z/(n(r, t), T(r, £)). Use of ( |3.10| ) then gives finally 

3A(r, t)(2A(r, t) + i/(r, t)) 2 = i/(r, t)a 2 . (3.13) 



This shows that X(n(r,t),T(r,t)) is the same as ( |2.16|) for uniform shear flow, except that 



the density and temperature are replaced by their values for the general nonequilibrium 
state. With this result for X(r,t) the solution ( |3.11| ) is completely determined. 



Next, consider the solution to ( |3.8| ) for the contributions to first order in the spatial gra- 
dients. As shown in Appendix [FJ, the right side is a linear combination of the hydrodynamic 
gradients. Consequently, p 1 ' necessarily has the same form 
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f^\V,y a (r,t))=X n ^ + X T ^ + X Uk ^ (3.14) 
where the coefficients Xp^(y a (r,t),V) are functions of the velocity to be determined from 



substitution of (|3.14|) into ( pT^ 



fd {0) , \ 

I - L(v ,a) + \6u- V V ' J X(3 >k - aX Uxjk 6(3 ;Uy = -Yp, k . (3.15) 

Here (3 denotes n, T, u x , u y , and u z , and i, k represents Cartesian coordinates. The last term 
on the left side originates from = ^r^r^ u k and Q3.9|) . The functions Yp^ are given 

explicitly in terms of and the flux t ■] , which must be determined self- consistently from 



The procedure is to solve (|3.15[) to find in terms of and then to use this result 
to derived a self-consistent equation for t\j . Further details can be found in Appendix 
along with the explicit solution to ( p. 151 ). 

In this way the complete normal solution near the uniform shear state is constructed to 
order e. The heat and momentum fluxes then are calculated from ( |3.4j ). The first terms, 
t\j\a) and q\ (a), represent the transport properties of the local reference state. They are 
the same as those of the previous section, fl2.18| ), except that the density and temperature 
are replaced by their values for the general nonequilibrium state. An important consequence 
of this replacement is that their gradients are no longer zero and thus they now contribute 
to the hydrodynamic equations. The second terms, 4j (a) and qp (a) provide new transport 
coefficients describing dissipation due to spatial deviations uniform shear flow 

It is understood that u, A(a), p, and the new transport coefficients all are functions of the 
local nonequilibrium temperature and density, as well as the shear rate a. The method 
for calculating 7^ (a)and C a ,j{ a ) anc ^ some detailed forms of the coefficients are given in 
Appendix |C[ a more complete listing can be obtained on request from the authors. Since 
the reference state is anisotropic there are new transport coefficients, reflecting the broken 
fluid symmetry, which do not exist for hydrodynamics near equilibrium. Fourier's law for 
the heat flux is modified by a thermal conductivity tensor, ^j(a), which has the form 

11 



characterized by five scalar "thermal conductivities" characterizing this tensor. The 
anisotropy is in the x, y plane so, for example, a temperature gradient in the x— direction 
leads to a heat flux in both x— and y— directions. At zero shear rate only the contribution 
from £<r (a) survives with £j^(0) = k, where k = 5fcep/2mz/ is the thermal conductivity 
for the BGK model. An additional difference from Fourier's law is a contribution to the 
heat flux from a density gradient characterized by the transport tensor, C n j( a )- This has 
a representation like (|3.17| ) in terms of five additional scalar transport coefficients. The 
asymptotic behavior for small shear rates is Cnj( a ) ~^ ( a ij + a ji) m {fl / v) 2 which vanishes 
for zero shear rate as expected. The momentum flux is expressed in terms of a viscosity 
tensor, 7^ (a), of degree four which is symmetric and traceless in ij. There are 19 inde- 
pendent viscosity coefficients. At zero shear rate %{{a) reduces to Newton's viscosity law, 
lirn( a = 0) = T)(5uSj m + S im 5j£ — ^SijSim), where r\ = p/u is the shear viscosity for the 
BGK-model. 

The corresponding non-linear hydrodynamic equations are obtained directly from these 
results and the conservation laws 

D t n + nV ■ 5u = 0, (3.18) 

-nknDtT - a-V — + —\t {0) - ^ ^1 
2 nkBVtI ailllm dr m + dr 3 ^ llm dr m J 

d dn dT 

-57^57 + ^57! = °- (3 ' 19) 

D,Su„ + P' l -^r\t l S - 7£^] + akM = 0, (3.20) 

where D t = d t +u • V. This is a primary result of our analysis. These equations are analogous 
to the non-linear Navier-Stokes equations for a fluid near local equilibrium, and reduce to 
them for zero shear rate. More generally, the reference state is local shear flow which can 
be very far from equilibrium. Furthermore, there is no restriction on the deviations from 
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this reference state, 5y a , since the non-linear dependence of all coefficients on y a (r,t) has 
been retained. The spatial gradients relative to shear flow must be small, however, and the 
equations are accurate to second order in these gradients. The terms proportional to 7^ and 
t\f in ( |3.19| ) represent viscous heating due to the excess gradient of the flow velocity relative 



to the reference state. There is no viscous heating from the reference state alone, since the 
local thermostat has been chosen to cancel it, but for states near uniform shear flow the 
thermostat cannot compensate for effects due to the gradients relative to the reference state. 
The implications of this new hydrodynamic description far from equilibrium are elaborated 
in the next section by an analysis of the associated linear hydrodynamic modes. 

IV. HYDRODYNAMIC MODES 

The above hydrodynamic equations are restricted to small spatial gradients relative to 
the reference state of uniform shear flow. If in addition the initial perturbations 5y a (0) are 
small then these equations can be linearized with respect to 5y a (t). The resulting set of five 
linear equations defines the hydrodynamic modes, or linear response excitations to small 
perturbations. If all of these modes decay in time the state is linearly stable. Otherwise, 
a growth of these modes signals an onset of instability that is ultimately controlled by the 
dominance of nonlinear terms. In this section we determine the hydrodynamic modes for 
states far from equilibrium and contrast them with those for states near equilibrium. 

The linearized hydrodynamic equations follow directly from (|3.18|) - ( |3.20| ) 



' d_ _d 

dt 13 3 dr; 



+ CLijTj——- ) 5n + n s V -5x1 = 0, (4.1) 



- + a l3 r 3 - \ 5T + — (-a lij>s + t^) — - — 



d 2 5n d 2 5T 



,n3,s dridrj Tj,s dridr 3 j 



'd d\_ _ x 

- + a ij r j — 8u k + Ps 



(4.2) 

(°)\ /«+(°)\ AXT 1 ^2; 



dt? k >\ d6n | / dtr k > \ d6T_ lk dH Ul 



dn dr< \ dT dr, nm ' s ' dr.dr 



+ a k j5uj = 0, 
(4.3) 
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To analyze these equations it is convenient to transform to the local Lagrangian frame, 
r[ — Ti — u si (r)t = Aij(t)rj. The Lees-Edwards boundary conditions then become simple 
periodic boundary conditions in the variable r'. A Fourier representation is defined by 



(M) = J dr'^'5y a (r,t) = J dre^Sy^t), 



(4.4) 



where the periodicity requires ki = 2riiii/Li, where n« are integers and Li are the linear 
dimensions of the system. In the second equality k(£), is given by 



ki(t) = kjkjiit). 



(4.5) 



The linearized hydrodynamic equations in this Fourier representation are, 
d 

— 5y a + (A au - ikj(t)B aU;j + k,j(t)ki(t)D auS )5y u = 0, 



(4.6) 



where in addition the dependent variables now have been scaled to dimensionless forms, 

-u > . (4.7) 



m 



5n 36T _ 

V V 27;' 



A summation convention applies and Latin indices denote Cartesian coordinates. The three 
matrices A a p, B a p, D a p are 



(4.8) 



B, 



af3j 



'ksTs 
m 








Ps \ dn 



n 

Ps 



dn 

'd£ 
dn 








in 

3 Ps 

lis. 

3 p s 

lis. 

3 p s 



dt 



(0) 



<9T 



dt 



(0) 



dt 



(0) 



dT 



Ojx 


3jy 




B2xj 




B 2z j 






























where B 2ij = \f\^;{-a^j,s + 4>)- The matrix D a/3 is 



(4.9) 
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D 



aPjl 



3 p s ^nl,B 3 p s ^ Tl , s 


















V 















— 1 ?2 ; —1 1'£ —1 72; 

P Txl,s P Tyl,, P %l,s 

—ljy — 1 jy — 1 jy 

P 1x1,8 P Tyl,s P %l,s 

— 1 jz — 1 72 —l 72 

p %i,s p r y i, s p %i,sj 



(4.10) 



The homogeneous solution to equations (|4.6|) can be calculated easily by putting k = 0, 



.(*) 



-At 



all 



5yp(0) = [l-M] aP Sy p (0). 



(4.11) 



The second equality follows from the property A 2 = 0. Consequently, all fields are constant 
except 5u x which behaves as, 



Su x (t) = u x (0) — at5u y (0). 



(4.12) 



The homogeneous state is unstable to an initial perturbation in 5u y , leading to unbounded 
linear change in time. Stability is still possible at finite k if this behavior is modulated by 
exponential hydrodynamic damping factors ~ e -^k 2 t w ^ a > q 

To simplify the analysis at k 7^ we allow the perturbation only along the velocity 
gradient direction, i.e. k =k y. In this case the linear hydrodynamic equations have time- 
independent coefficients (i.e., k(t) = k), 



Sy v = 0, 



a/3 



A af3 - ik B a/3 + k 2 D af3 , 



and the matrices B and D now take the simpler forms, 



B, 



a/3 



k B T s 
m 







at 



(0) 



Ps \ dn 






2T S ( dt 
3 P. 



E>2xy B2yy 



dt 



(0) 



dn 



2 J-s I UL VV 



3 p a 



9T 














(4.13) 



(4.14) 
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(4.15) 



zy,s J 



Equation (gig) can be solved by Laplace transformation, 



(4.16) 



6y a (k,z)= / rfte- t ^ Q (k,t) = [^ + F(k,z)]-^(k,t = 0). 

The eigenvalues u^(k,a) of the matrix F(k,z) define the five simple hydrodynamic poles 
at z = — o;W(k, a), which determine the dominant dynamics of the §y a (k,t) at large t and 
small fc. At equilibrium (a = 0), the hydrodynamic modes of the Navier-Stokes equations 
are recovered (two sound modes, a heat mode and a two-fold degenerate shear mode) for 
long wavelengths (k — > 0), 

\ 



V 



icfc + r/c 2 

D T k 2 
(v/p) k 2 
(V/P) k 2 



(4.17) 



/ 



where c 



5/3/3m is the sound velocity, T = Z>r/3 + — 



— 1/Pmu is the sound damping 
constant, D T = £,Ty,s( a = 0)/pCp = l//3mv is the heat diffusion coefficient, C p is the specific 
heat per unit mass, and rj/p = l//3mu is the kinematic viscosity. The equilvalence of T, 
Dt, and rj/p is a peculiarity of the BGK model. These coefficients are positive so that 
Eqs.( |4.17| ) represent damped excitations. Corrections to these dispersion relations are of 
order k 3 , describing an expansion that is analytic in k about k = 0. 

For finite shear rate, the modes are more complicated and the behavior at long wave- 
lengths is qualitatively different. To be more precise, consider the case of k — > at fixed, 
finite a. It follows directly from ( [4 . 1 6| ) that the hydrodynamic modes have the asymptotic 
behavior 
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c x [a)k 2 

+ i v/3)c 2 (a)A: 2 / 3 + |(1 - i \/3)c 3 (a)A; 4 / 3 + c 4 (a)A; 2 
-|(1 - ^ v / 3)c2(a)A; 2//3 + ±(1 + i \/3)c 3 (a)/c 4 / 3 + c 4 (a)A; 2 
c 2 (a)A; 2 / 3 + C3(a)A; 4 / 3 + C4(a)A; 2 
(r;(a)/p)fc 2 



with the coefficients Cj(a) given by 

'9n(a)dm{a) d V (a)dt^(a) 



ci(a) 



3T <9n 



<9n <9T 



«9T 



(4.18) 



(4.19) 



c 2 (a) 



2a2 M + «w) Sf » ,(a) 



3n 2 mkB 



dT 



1/3 



(4.20) 



c 3 (o) 



9n 2 mkBC2(a) 
1 



3nmc 2 (a) 



S^'-a^a) + 4?(a)) - a (77(a) + ^(a))^* 



dT 



<9T 



(4.21) 



04(a) 



ci(a) 



p- 1 fe(«)+7») + 3^^(a 



(4.22) 



The modes (|4.18|) represent two oscillating modes and three purely damped modes, just as 



in the Navier-Stokes case (|4.17|) . However, there are two important qualitative differences. 
First, the long wavelength behavior is non- analytic in k about k = and is given by a 
power series in k 2 ^ 3 . Thus, for example the purely damped modes do not represent diffusive 
behavior as in the Navier-Stokes case. This non-analytic behavior with respect to k is due 
to the fact that the reference matrix at k = is not diagonalizable and the eigenvalues are 
not analytic about k — a — 0. Therefore, recovery of the form of the modes near equilibrium 
requires that k and a be taken to zero in a related way (see below). The hydrodynamics 
for an alternative choice of thermostat, discussed in Appendix A, has dispersion relations 
that are analytic about k = a = 0. A second critical difference between( [4.18] ) and ( [4.17] ) 
is that the two oscillating modes are unstable in the long wavelength limit because the 
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coefficient 02(a) is positive for all a. This means the modes include excitations that grow 
in time. Eventually, the deviations 5y a grow beyond the limitations of the linear equations 
and the full non-linear hydrodynamic equations are required to determine their ultimate 
values. These will be different from those of the reference state, representing the fact that 
the reference state itself is unstable. 

It is possible that the hydrodynamic modes are stable at shorter wavelengths. This is 
in fact the case, as can be seen by solving Re[o;W(k s , a)} = to determine the stability line 
k s (a) in the k-a plane. This is illustrated in Figure where the calculation was performed 



using the exact eigenvalues rather than the small k expansion of (}4.18|) . Dimensionless 
variables are used, k* = kl and a* = ar, where r = v^ 1 is the mean free time, i = Vq/u is 



the mean free path, and Vq = y2kBT/m is the thermal velocity. Above this line the modes 
are stable, while below this line they are unstable. This prediction of a long wavelength 
instability has been verified quantitatively by comparison with Monte Carlo simulation of 



the same kinetic equation from which this hydrodynamics was derived |TT| . Further analysis 
of this instability and comparisons to simulations is reported in detail elsewhere |12|| . In 
the following we focus on the stable domain of Fig. [TJ. To study the stable dynamics a new 
dimensionless variable, x = k*/a*, is introduced. The hydrodynamic modes are expressed as 
functions of k and x, i.e., tu^(k, a) = o;W(k, x), and the expansion about k = is performed 
at fixed x. Physically, this involves controlling both the shear rate and the wavelength 
simultaneously to assure that the system is stable (sufficiently large x) while approaching 
the long wavelength limit. To simplify the calculation, a system of Maxwell molecules is 
considered (interatomic potential ~ r -4 ). In this case z/(n, T) and X(n,T) are independent 
of the temperature. The hydrodynamic modes for the stable domain are then obtained for 
x > ;4r, which lies above the dashed line shown in Figure [I] 
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1 ick + Tk 2 -la 2 ^ 
-ick + Tk 2 - fa 2 



10 



D T k 2 



B T k 2 + fa 2 



(4.23) 



y (V/P) k 2 j 

For this expansion at fixed x the eigenvalues are again analytic functions of k and can be 
interpreted as perturbations of the Navier-Stokes modes due to small but finite shear rate. 
The restrictions on x imply < k* < 1, so the shear rate dependence is small but not 
necessarily in the Navier-Stokes limit. More generally, the entire stable domain including 
larger shear rates can be accessed for x s = k/k s (a) > 1 and evaluating cZ;(k, function 
of k exactly. 

The linearized hydrodynamic variables, 5y a (k, t), can be expressed in terms of the eigen- 
values and eigenfunctions, 



a (k,t) = ^e- W ( k ^ ) (k,x)(^)t(k,o;)%(k,0) 



(4.24) 



where {C^} are the eigenvectors and {rj^} are the associated bi-orthogonal set defined by 
^aiVa^YCa ^ = To illustrate the effects of the shear rate the analytic results for small k, 
( [4.23j ), will be used. The corresponding eigenvectors are 



( 



3 2 i n 

5' V5' V5i' ' 



™ 10 

5' V 5' V5i' ' 



§, 1,0, 0,0 



V V 3' X ' V&n' U ' U 

(0, 0, 0, 0,1) 



\ 



(4.25) 
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V2 V 5' TlO' ' 2' u 



(- 



0. 



0,0 



3/3 3 
.5V 2' 5 



2' ' 2 ' ' 
i\/3a 



0,0 



V 



2 ' "> 

(0, 0, 0, 0,1 



(4.26) 



where it is understood that k is restricted to the stable domain. Then the response of the 
density to an initial density perturbation is found to be 



5n(k,t) 



- e -( rfc2 -lr)*cos(A;ct) + e" 



D T ri 



3 4g 2 1 

—e 5" 
5 



Sn(k,0) 



(4.27) 



Thus, for small shear rates there is an enhancement of the amplitudes for both the sound 
and heat modes. All other hydrodynamic variables can be calculated in a similar way from 
( |4.24j ) as well for both the stable and unstable regions, including small and large shear rates. 



V. DISCUSSION 

The objective here has been to study transport far from equilibrium for the special 
nonequilibrium states near uniform shear flow. Dynamical properties of states far from 
equilibrium are not well understood due to their complexity and technical difficulties with 
the formal theories of nonequilibrium statistical mechanics. However at low density the 
Boltzmann kinetic theory provides a controlled formulation of this problem. There are still 
difficulties for practical applications so a kinetic model has been used to allow a detailed 
analysis for a special class of states near uniform shear flow. First, the BGK-Boltzmann 
equation was solved exactly for the steady state distribution at uniform shear flow, and the 
corresponding transport properties given as a function of the shear rate. Next, a solution to 
the kinetic equation was obtained for a class of states deviating from uniform shear flow by 
small spatial gradients in the hydrodynamic fields, using a variant of the Chapman-Enskog 
approximation method. By "variant" is meant that a local form of the stationary solution 
for shear flow is used as a reference function rather than the local equilibrium distribution 
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function. In general this reference state is very different from a Maxwellian and can be very 
far from equilibrium if the shear rate is large. The Chapman-Enskog expansion was used to 
determine the distribution function to first order in the gradients. The result is a "normal 
solution" for which all space and time dependence occurs through the hydrodynamic fields. 
These fields must be determined from hydrodynamic equations which follow from the exact 
conservation laws. The irreversible parts of the hydrodynamic fluxes were determined as 
functions of the hydrodynamic fields and their gradients using the normal solution with 
results of the form £*■ = ~V^f ~ i& for the momentum flux and q* = - C n ,j fr 

for the heat flux. Since these coefficients are calculated near the stationary state of broken 
symmetry there are many new transport coefficients (7^, £y -, •), in comparison to the case 
of states near equilibrium, which depend on the shear rate. With these new expressions for 
the fluxes, a closed set of the generalized hydrodynamic equations was derived near uniform 
shear flow. In summary, a complete description at both the kinetic and hydrodynamic levels 
has been given for a wide class of states arbitrarily far from equilibrium. In particular, 
the analysis provides a rare example of the relevance of a hydrodynamic description far 
from equilibrium. The corresponding hydrodynamic modes were calculated to order k 2 for 
arbitrary shear rate from the linearized hydrodynamic equations. An unexpected result is 
an instability for any finite value of the shear rate at sufficiently long wavelengths. This 
prediction has been verified quantitatively by comparison with direct Monte Carlo computer 
simulation of the kinetic equation, confirming the validity of the hydrodynamic description 

mil- 



The hydrodynamic analysis was carried out here only for spatial variations along the 
velocity gradient. More complex dynamics is expected in the general case of arbitrary 
direction for the spatial perturbation and a description will be given elsewhere. In addition 
to the stable and unstable exponential time dependence of the hydrodynamic modes, there 
will be algebraic modulations due to the fact that the hydrodynamic matrix at k = cannot 
be diagonalized. This can lead to initial growth of perturbations of uniform shear flow even 
when the hydrodynamic modes are stable. 
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The approach taken here can be extended to several other physically interesting ref- 
erence non-equilibrium states such as a constant temperature gradient, or combined heat 
and momentum transport. Exact solutions to the BGK kinetic equation are known for these 
cases H so the reference distribution for the Chapman- Enskog expansion is available. Other 
directions of extension include applications at higher densities. Recently, a BGK-like kinetic 
model for the dense fluid Enskog kinetic equation has been described and applied with suc- 



cess to shear flow ||22|| . This provides theoretical access to densities relevant for molecular 
dynamics simulations. 

Although attention has been focused on transport the analysis can be extended in a 
straightforward way to describe fluctuations in uniform shear flow. The reason for this is 



the close relationship of the kinetic equations for fluctuations to that for transport |23| , |24 
For example, the kinetic equation for phase space fluctuations at two times is governed 
by the linearization of the kinetic equation for transport. As a consequence, the linear 
hydrodynamic equations studied in Section 4 also can be used to compute the hydrodynamic 
part of time correlation functions such as the dynamic structure factor measured in light 
scattering. Recently, kinetic models for fluctuations have been developed that are self- 
consistent with the BGK kinetic model for transport [2~5|. Their detailed application to 



fluctuations in uniform shear flow, including anomalous long range spatial correlations |26| 
will be discussed elsewhere. 
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APPENDIX A: THERMOSTATS 



To establish a steady state for uniform shear flow a non- conservative external force is 
introduced in section |I| to compensate for viscous heating generated by the Lees-Edwards 
boundary conditions. In section [TTT] for states near uniform shear flow, a local form for this 
non-conservative force is used 

F ext (r, t) = -mX(n(r, t), T(r, t))(v - u(r, t)). (Al) 

With this local form it is possible to impose df^T = 0, which led to the results that 
A(n(r, t), T(r, t)) is the same as A s (n s , T s ) at uniform shear flow given by (|2.17f) , except that 
the temperature and density are replaced by those for the general nonequilibrium state. The 
resulting hydrodynamic equations express viscous heating only due to the gradients relative 
to uniform shear flow. An alternative choice is to use the simpler case of the constant 
X(n s ,T s ) in ( |A1|) even for states near uniform shear flow. The advantage of this choice is 
a simpler implementation of the Monte Carlo simulation method for the solution to the 
kinetic model equations. Obviously, it is consistent with the conditions for the steady state 
uniform shear flow. However, this choice only partially compensates for the viscous heating 
at local uniform shear flow, i.e., df^T ^ 0. In this appendix, the changes from the results 
of section [HI] in both the Chapman-Enskog solution and the hydrodynamic equations due 
to this alternative choice are described. 

To lowest order in the uniformity parameter the BGK equation with the new thermostat 
becomes, 

[ — - L(v', a) + A s 6u • V V , /(°) = uU (A2) 

Here and below it is understood that A s = A(n s , T s ) is independent of space and time. Since 
the solution is normal, its time dependence occurs only through the hydrodynamic variables 
and the contributions from the time derivative, , can be calculated using the conservation 
laws to lowest order in the uniformity parameter, 
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= 0, + a tJ 6 Uj (r, t) = 0, (A3) 

^n(r,t)k B - )T Jf ,t \ a lJ t\f(r,t) + 3\ s n(r,t)k B T(r,t) = 0. (A4) 

These differ from the results of the section [TIT] because it is no longer possible to choose X s 
to compensates for the viscous heating, i.e., Oijt\j\r,t) ^ — 3 A s n(r, t)k B T(r, t). Instead, 
Eq.(IA3) becomes, 



■(oijtif + 3X S nk B T)-^= - L{{V - 8u),a) 



f {0) = vfe, (A5) 



3nk B " lJ 13 ' " OT 

where the temperature derivative term is new. To solve the Eq. (|A5| ) , we need to know 
the temperature dependence of the zeroth momentum flux, tj®. For power law potentials 
this follows from the temperature dependence of v. The analysis is simplest for Maxwell 
molecules where v is independent of the temperature. Then tj® ~ T (as is verified a 
posteriori) and the solution to Eq.( |A5| ) is 

( 771 \ ^ poo 

f^\V,y a (v,t)) = u(r,t)n(r,t) [————} / dre^^e 3 ^ 



27rk B T(r,t) J Jo 

exp (- 2 ^y (r;t) e2w(r ' t)T [ A ^(- r )K- - *))]') - (A6) 

w(n(r, t),T(r, t)) = at £fc*) (A7) 

V V ' ; ' K ' JJ 3n(r,t)k B T(r,t) V ; 

The solution is still only implicit since t^(r, t) , or equivalently w(r, t), must be determined 
self-consistently using (|AB|). The result is, 



3w(r,t){2w{r,t) + z/(r,t)) 2 = z/(r,t)a 2 . (A8) 

This equation is the same as that for A(n(r, t), T(r, t)), Eq. ( |3.13| ), showing that w(r, t) has 
the same functional dependence on a, u(r,t), and T(r,t). This proves that the zeroth order 
solutions to the BGK-Boltzmann equation, ( |A6| ) and ( |3.11| ), for the two different thermostats 



are the same. The lowest order equations for the temperature differ, however. Use of ( [AT] ) 
in (A"4) gives, 
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dj®T(r,t) = 2(w-\ a )T(r,t), 



(A9) 



which vanishes only at the steady state. 

Now consider the first order solution in the Chapman-Enskog expansion, 



'9(0) 
dt 



L(v', a) + A s 5u ■ V v ' / 



(i) 



and look for solutions of the form, 



/ (1) (v', y a {v, t)) = X n>i (v', y a (r, t))^- + X T>i (v', y a {r, t))%- + X Uk ,(V , y a (r, t)) 



dn 



dT 



(A10) 



dr, 



35uk 



(All) 



The equations for the coefficients, X/3 t i(v';y a (r,t),), are determined from substitution of 
into the ( |A10| ) and they are 



dt 



— L(v' a) + A s 5u ■ V v ' 



X n ,k 
Xr,k 

X u - b 



+ 



2 (A, - w) X TM 
adi v X u ,b 



(y \ 

*n,k 



(A12) 



V Iu J> k j 

where i,j, k represent Cartesian coordinate. The functions Yg^ are same as those given 
in Appendix B. Since the equations for X n> k,XT,k are different from ( |3.15| ), the transport 
coefficients for the heat flux, £^ ■, £^ •, are also different for this new thermostat. However, 
the equation for X Ui ^ is unchanged so the transport coefficients for the momentum flux, 7^, 
are unchanged. The corresponding hydrodynamic equations are the same as in the section 
IIJ, except for the temperature equation, which becomes 



3 , „ „ u 05ui , d5ui (o) ij OSui . 



-nk B D t T - aij^ 



y 

m dr 



+ 



dr 3 1 13 



\t { ' - ^ 



dr r 



(A13) 



The new term, 3 (w — A s )nfcsT, on the right side represents the viscous heating due to 
spatially uniform deviations from the steady state. 
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APPENDIX B: CHAPMAN ENSKOG EXPANSION 



In this appendix, Eg. (|3.8|) for /W i s solved, 
f d (o) 



Of 

)(i) _ am 



L(v', a) + X5u ■ V^J /« = - (A (1) + V • V) /«», (Bl) 



where D t = + a^rj-^. The right side of Eq.([Bl[) can be written in terms of hydrody- 
namic derivatives, 

(A (1) + v ■ V)/(°) = £ ^(A W + v • V)y a , (B2) 

a=l 

where y a = {n(r, t), T(r,t), <5u(r, £)}. The time derivatives can be replaced by first order 
spatial gradients of hydrodynamic variables using the corresponding hydrodynamic equa- 
tions obtained from Eqs.( [2.6|) - ([2.8|) 

D ( pn + V • (ndu) = 0, (B3) 
\nk B {D? + Su • V)T + aij tf) + 4 0) ^ + V • q<°> = 0, (B4) 

0,(0) 

D^Sm + 5u ■ Vdui + p" 1 ^- = 0. (B5) 
Use of these equations in Eq.(|B2"D gives 

(A (1) + v ■ V)/<°> = K nii (v';y Q (r,0)^ + Y T ,i(v'; y a (r, + Y Uk;i (v';y a (v,t))^. 

(B6) 



The coefficients, Y^, are given by 

Y n>i = ——{Vi - 5ui) - p — — , (B7) 

an ooui on 

^-^ W -^)-p-g^ (B8) 
^ = -n^H* - (*2? + ^-^Irl + - *0, (B9) 
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where i, I, k, m denote x, y, z. The derivatives of f(°' can be obtained simply from Eq.( [3.11| ), 
d f{0) f(o) fdu\l 



dn n 



+ (s)>-^i + (s)^- (B10 » 



Here the functions /i Q are defined by, 

h u = j Q drre~ , ' T+3XT f £ (e XT A tJ (-r)(v' J -5u,)), ( B1 3) 



x/^A^-r)^-^)), (B14) 
^ = / dTe-^V^A^-T)^ - ^OlV^e^iC-r)^. - <^)), (B15) 

« 



CXD 

= J rfre-^ +3A V A ^ fc (v / ,r)/,(e A ^(-r)(^.-5^)), (B16) 

where the functions g%(v',t) in the expression for h Uk are 

= ( v x - Su x) + at(v' y - 5u y ), 
g y = at(v' x - 8u x ) + (1 + a 2 t 2 )(v' y - 5u y ), 

g z = v' z - 5u z . (B17) 

The right side of Eq. (|B1|) now can be written in terms of the coefficients, Y a ,i, and the 
hydrodynamic gradients. Note that the Y a ^ are given only implicitly since they depend on 
the unknown flux . The procedure is therefore to determine in terms of tjj , and then 
use that result to obtain a self-consistent equation for t\j\ 

Clearly, f^ 1 ' must have the same form as the right side of flB6|) 
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/« = Xn^V, y a (r, t))j£ + X T4 (V, y a (r, t))^- + X Uk ,(V , y a (v, t))^. (B18) 



Substitution of this form into flBl| ) gives the equations 

\-QjT - L{v , a) + A 5u i-Q-j ) x a,k - aX Uxtk 5 a>Uy = -Y a>k , (B19) 

where a = n,T, u x , u y , u z . These equations have a form similar to Eq.(4.8) and can be solved 
in a similar way to get 

roc 

X a , k = J o dre-^e^lad^X^ie^A^-r)^ - 8 Uj )) 

-Y^e^A^-r)^ - 8 Uj ))]. (B20) 

This implicit result is sufficient for determination of the transport coefficients, as described 
in the next section. 

APPENDIX C: TRANSPORT COEFFICIENTS 



In this appendix, the transport coefficients defined in Eqs.( [3.16D are calculated using the 



solution to Eq. (|3.15| ). First, the coefficients from momentum flux are considered, 

lira = ~ J dVmCiCjXu^m, (CI) 

where c = v' — 5u. As an example, 7^ is analyzed in detail. Substitution of the solution, 
X Ux>y , (|B20| ) into the definition of 7^, flC]} gives 

/roo 
dv'mc x c y J o dte-^Y^e^A^-t)^). (C2) 

To proceed it is necessary to know the temperature dependence of v and A. The simplest case 
is that of Maxwell molecules for which v and A are constants. Then 7^ can be calculated 
with a sequence of change of variables: first v' — > v' + Su, then v' — > e~ xt V, and finally 



v\ — > Aij(t)v'j. In this way 7^ is found to be, 
7^ = J™ dte- (u+2X)t m J dw (v x - atv y )v 



^V + a, 7-) + ^ 

\ l xy 1 a lm Jxy I < or u y 



3nk B 8T K xy ~" ni *v> 35u 
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Use of the result (|37T2|) for tg and the relation (|37T3|) for A gives the final value of the 7xy > 

r y = v{v - 2A) P . (C4) 

The same method can be applied for all the transport coefficients, 7^, • and £y •. A 
complete listing of the results can be obtained from the authors on request. 
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FIGURES 

FIG. 1. The solid line separates the stable (above) and unstable (below) domains in the k* = k£, 



= ajv parameter space. The dotted line is the asymptotic form obtained from section IV 
> 2/V5). 
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